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Abstract 

We develop a numerical method for realizing mean curvature motion of in- 
terfaces separating multiple phases, whose areas are preserved throughout 
time. The foundation of the method is a thresholding algorithm of the 
Bence-Merriman-Osher type. The original algorithm is reformulated in a 
vector setting, which allows for a natural inclusion of constraints, even in the 
multiphase case. Moreover, a new method for overcoming the inaccuracy of 
thresholding methods on non-adaptive grids is designed, since this inaccuracy 
becomes especially prominent in area-preserving motions. Formal analysis of 
the method and numerical tests are presented. 
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1. Introduction 

This work develops a method to compute the length-shortening evolution 
of interfaces between an arbitrary number of phases in arbitrary dimension 
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under the constraint that the area of each phase is preserved throughout time. 
Such evolutions often appear in situations where interfaces move according 
to their geometry, while the mass of each phase remains constant (e.g., grain 
boundaries in ternary alloys, crystal growth, multiphase flows or formation of 
soap film bubbles). In these examples, the motion is driven by the decreasing 
energy of the internal interfaces, which are out of equilibrium. This kind of 
motion also has applications in image processing (denoising, segmentation), 
in biology (modelling of vesicles and blood cells), in the description of isolated 
gravitating systems in general relativity ^| , and other research fields. 



Strictly speaking, since area preservation is a global constraint, one can- 
not consider a constrained curvature flow directly. Therefore, we have to 
start from a more basic aspect of the motion, such as the energy. In par- 
ticular, we consider the constrained steepest descent of the "length energy" 
of each interface, which counts the measure of interfaces weighted by their 
corresponding interfacial tensions. The steepest descent of the length energy 
without any constraint gives the classical mean curvature flow. On the other 
hand, in the case of two phases, the area-constrained gradient flow of this 
energy corresponds to evolution by mean curvature, minus a time-dependent 
term (equal to the average mean curvature over the interface). The situ- 
ation is analogous for more than two phases but the nonlocal term has a 
complicated form which depends on the configuration of each interface. 

The subject is also mathematically interesting, because it is one of the 
most simple problems with nontrivial limiting behavior. It is well known 
that mean curvature flow shrinks uniformly convex smooth hypersurfaces 
smoothly to a point in finite time. On the other hand, the area-preserving 
mean curvature flow converges to the solution of the isoperimetric problem, 
i.e., a sphere j2|, S, 0, Isf . However, the area-preserving flow may drive general 
embedded hypersurfaces to self- intersections, as was shown in On the 
other hand, [3| and jsj proved that if the initial surface is sufficiently close 
to a sphere then it converges to the sphere even if it is not initially convex. 
Due to the complexity of the multiphase case, there are only a few results 
concerning the stability of junctions under area-preserving flow, see j^, [l^ 
and the references therein. 

Since evolution of surfaces is an intensely studied subject of practical 
interest, a number of analytical and numerical methods have been developed 
to treat motion by mean curvature. Many of these methods can be applied to 
the constrained motion addressed here; let us summarize the known results 
with emphasis on the multiphase case and volume preservation. 
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Perhaps the most basic approach is to use the definition of the motion 
itself. That is, in the two-phase case, one computes the evolution of the 
interface directly from its velocity: 

v{x) = {—k{x) + Ka)n{x), a.e. x G dP{t), 

where P{t) denotes the region occupied by one phase, k is the mean curvature 
and is the average mean curvature over the whole interface. Algorithms 



based on this idea are called front tracking methods [111 . |12| | . They directly 
approximate the interface based on the Huygens' principle and are effective 
for computing the evolution of smooth surfaces without topological changes. 
Although this method is simple in principle, if interaction of different parts of 
the interface occurs, a complicated decision algorithm is necessary to proceed 
with the computation, and this becomes increasingly more involved in higher 
dimensions. Higher dimensions and a larger number of phases also complicate 
the calculation of curvatures and their averages over the interfaces. 

A more general framework is provided by the level set approach which, 
thanks to its implicit representation of the interface, is able to deal with 
topological singularities and nonsmooth data. In this approach, the initial 
interface 5P(0) is expressed as the 0- level set of a function 0(a;, 0), and the 
mean curvature flow is achieved as 

dP{t) = = 0}, 

where (j) is the viscosity solution of the Hamilton- Jacobi equation 

The constrained flow can be realized in this setting by considering the area- 
constrained gradient flow of the length energy functional written in terms of 
the level set function 

L{dP{t)) = j 5{(t){x,t))\V(t){x,t)\dx, 

where 5 denotes the Dirac delta function. It is still necessary to calculate the 
curvature values, but this method can be extended to the multiphase setting 
by introducing as many level set functions as there are phases and imposing 
an additional constraint so that the level sets do not overlap or create vacu- 
ums (see 13| or However, such a constraint has an unwanted impact 
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on the flow |15l] and the phase areas are not adequately preserved during 



the computations. The first problem was solved in [15| and [16| by employ- 
ing signed distance functions. Area-preserving motions are also addressed in 
these works but are limited to two-phase flows. Another multiphase mod- 
iflcation of the level-set method, which has some similarities to our own 



approach, was developed in [1^] , where the constraint of [13[ was replaced by 
a projection step. However, the impact of the projection on the dynamics of 
the interface was not analyzed. 

The area-preserving mean curvature flow arises as a limit of the following 
nonlocal mass-preserving diffusion equation 18. 19. 2ol 



Ut = Au- \w'{u) + [ W'{u) dx, 

where W is a double-well potential and e is a small parameter related to 
the width of the diffuse interface. It has been shown that, under suitable 
conditions, the set 

Peit) = {x; u,ix,t) > i} 

approximates P{t) with error 0(£:^| log5p). 

Based on this fact, the so-called phase fleld methods represent interfaces 
by thin layers in the solution and thus the resolution of this internal layer 
requires a very flne mesh. On the other hand, this approach handles topo- 
logical changes without trouble and does not require explicit computation of 
curvatures. An interesting computational analysis for the multiphase case is 
given in 
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The basic idea of this paper is to use another approach, often referred to as 
a thresholding method. We adapt the so-called BMO algorithm (sometimes 
also called the MBO algorithm), which was discovered by Merriman, Bence 



and Osher in [22|, to generate multiphase area-preserving motion. As far as 
is known to the authors, the existing works (except do not treat the 
approximation of multiphase length-shortening flow under area constraints. 

The BMO algorithm exploits the fact that short-time diffusion of the 
characteristic function of a region enclosed by an interface (i.e., its convolu- 
tion with the Gaussian kernel), evolves the interface according to its mean 
curvature. More precisely, the characteristic function of a region is evolved 
for a short time by the heat equation and then a thresholding step is carried 
out to obtain the new interface (given by the 1/2-level set of the diffused 
function). The main advantage of this approach is that it naturally treats 
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topological changes, produces no intercalary regions and does not require 
explicit computation of curvatures. Moreover, it is numerically attractive 
because of its stability and low computational complexity. 



This thresholding method was applied to multiphase flow in [22|, while 



23[ uses the BMO algorithm for constructing two-phase area-preserving cur- 



vature flow. However, the latter method lacks theoretical support, and it is 
not clear how to extend it to more than two phases and to more general mo- 
tions. We therefore introduce a different approach, which is also based on the 
BMO method. Our method can treat any number of phases in any dimen- 
sion and can be extended to more general motions, such as mean curvature 
motion with transport. 

The difficulty of the multiphase constrained flow is that the phases in- 
fluence each other not only locally via the shape of their interface, but also 
globally via their areas. The idea used to overcome this complication is to 
formulate the original multiphase BMO method in a vector-valued fashion 
and to realize the area constraint by considering a constrained gradient flow. 
This constrained flow presents a computational difficulty, due to the fact 
that the interfacial velocities are slower when compared to the flow without 
area preservation. That is, since the interface must move at least the dis- 
tance of the grid size at each time step, this places unreasonable restrictions 
on the grid resolution used in the numerical implementation. We are able 
to overcome these restrictions by introducing a technique of temporary and 
localized reflnement. 

The paper is organized in the following way. In section 2, we introduce 
the area-preserving mean curvature flow, as well as the BMO algorithm for its 
approximation. We discuss the numerical algorithm in section 3, and section 
4 concerns its implementation. Section 5 presents a number of numerical ex- 
amples and analyses of errors and model parameters. The appendix includes 
a number of theoretical results requiring technical computations. 



2. BMO algorithm for area-preserving mean curvature flow 

2.1. Area-preserving mean curvature flow 
2.1.1. Two-phase case 

We flrst fix the meaning of the terms "area" and "length" used in the 
sequel. Working in m-dimensional space, the word "area" shall mean the 
m-dimensional Lebesgue measure of the region corresponding to a phase, 
i.e., in 3 dimensions it regards the volume. The word "length", on the other 
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hand, shall always refer to the (m — 1) -dimensional Lebesgue measure of 
the boundary of a phase region, thus in 3 dimensions corresponding to the 
surface area of the interface. 




n 



Figure 1: Two phases divided by an interface. 



Mean curvature flow is related to systems whose energy depends on the 

length of their surface, e.g., soap films. We will explain the basic equations 
for curves in the two-dimensional plane, since the derivations in higher di- 
mensions are similar but lack transparency. Accordingly, let us consider a 
smooth Jordan curve 7 contained in a subregion as in flgure 1. Let the 
curve be parametrized: 



and in such a way that the enclosed region P is on the left side of the curve. 
Then the length of the curve is given by 



The gradient flow of the above energy can be found from its flrst variation. 
That is, for any smooth closed curve (p{s), s e [a, b] we compute 



where k is the curvature, n is the unit outer normal to the curve 7 at a given 
point and where we integrate with respect to arc length I: 



7(s) = (7i(s),72(s)), se[a,b], 



7(«) = lib) 





(1) 
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In the general m- dimensional case one obtains the same formula as ([T]), where 
K means the trace of second fundamental form divided by m. Hence, the 
fastest shortening of the curve occurs for the flow with normal velocity equal 
to minus (a multiple of) the mean curvature. 

Next we consider the curve-shortening flow under the constraint of area 
preservation. The area functional for region P reads 



and its first variation is 

A. 

de 



^ j{x,y)-ndl = ^ j (7i72 - 7271) ds 



A{-i + eLp)\e=o = / n-ipdl. 



Analogously, in any dimension the variation turns out to be the normal vector 
at each point of the hypersurface. 

Following the construction in 2J], we introduce a time- dependent La- 
grange multiplier A(t) for the area constraint and express the velocity of the 
area-constrained mean curvature flow by 

V = (—K + \)n. 

A precise expression for the the multiplier A can be obtained in the following 
standard way. From the fact that the area is preserved one has 



-A(7(t))= / v{t)-n{t)dl = 0. 
at J u) 



hit) 

Hence it follows that the integral of — k -|- A over the curve 7 vanishes at each 
time. This yields 

m = jA^ [ <t) di. 

That is, the Lagrange multiplier expresses the average mean curvature along 
the interface. 



2.1.2. Multiphase case 

We briefly derive the velocity of interfaces moving by area-preserving 
mean curvature flow in the multiphase two-dimensional setting. We assume 
that the number of phases k is finite and that the interfaces between different 
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pairs of phases form a finite collection of arcs. The boundary 7^ of phase 
region Pi can then be written as the union of the interfaces from all other 
phase regions: 

Here 7jj denotes the interface between phases Pj and Pj, and the index I 
expresses the fact that each interface may have several disjoint parts. In the 
following, however, we omit the decomposition using index / since it has no 
influence on the computations. 

Each interface 7'^ is considered as an oriented curve, which has the region 
Pj on its left side. For any point interior to the interface 7ij we define the 
normal n as the unit vector pointing into the phase with larger index, i.e., 
n is the outer normal to Pj if i < j. We remark that some of the curves 7jj 
may be empty. 

The energy of the system, considering arbitrary surface tension Tij for 7^- , 

is 

k 

dl- (2) 

i=l j>i J Hi 

This value is to be minimized under the condition of constant areas, that is. 




= const, i = 1, . . . , A; — 1. 



Introducing Lagrange multipliers A,, i = 1, . . . , k — 1, for each of the con- 
straints, the constrained variation reads 

i=l j>i ''Hi 1=1 \j>i ''Hi j<i ''"la J 

where v?*-^, i, j = 1, . . . , /c, denotes a smooth perturbation within 7jj. 

From this it follows that the magnitude of normal velocity v'^^ for interface 
lij {i < j) in the direction of n will be 

V"-^ = -TijK, + A, - (1 - 5jk)Xj. 

Here 8jk denotes the Kronecker delta, which arises from the redundance of 
the area constraint for phase P^ (it follows from the constraints on the other 
phases and the fixation of the domain). 
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Explicit representations of the Lagrange multipliers can be obtained as 
follows. Let us denote by L.^ the length of interface jij, by Li the length 
of the boundary of phase region Pi, and by Kij the average mean curvature 
through interface jij times the tension Tij: 

Then the preservation of phase area Pj gives the condition 



= J2 [ v'Ul-J2 [ 

j>i '•''J j<i Tj» 

^ / {-Tijn) dl-^ {-TijK) dl 

j>i "^Ty j<i ''iji 

j>i '^'Tij j<i "'T: 

+ J2(^-Sjk){-Xj) [ dl-Y,\ I 



j>i f'J j<i ''Iji 

dl 

j>i ''^ij j<i ''^ji 

— ~ ^ LijKij + ^ LijKij + Aj ^ Lij — — Sjk)XjLij. 

j>i j<i j^i j^i 

Since the above holds for i = 1, . . . ,k — 1, we obtain a system of linear 
equations for Ai, . . . , Ajt_i: 

fc-i 

LiXi — ^ ^ LijXj — ^ ^ LijKij — ^ ^ LijKij, i — 1, . . . , k — 1. 

j=l j>i j<i 

For a given configuration, the solution to this system gives the Lagrange 
multipliers. For example, in the case of three phases, as in figure 2, one has 
the velocities: 

v'' = -ri3«+(l-^)-r^-(l-^)-^% 

Here a — L13L23 + L23L12 + L13L12 and k^*, i — 1,2, 3, represent the average 
mean curvatures along the whole boundary of each phase, weighted by the 
surface tension. 
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Figure 2: An example of a three-phase configuration. 



When phases Pi and P2 are separated, i.e., when L12 = 0, the above 
formulae reduce to the form of the two-phase flow from section 12.1.11 The 
extension of the above calculations to higher space dimensions is natural, as 
we use only the notions of interfaces and their oriented normals. 

2.2. BMO algorithm 

For the sake of clarity, we explain our method in three successive steps. 
First we summarize the original idea of the BMO algorithm for two phases 
when no area constraint is present. Then we describe the existing algorithm 
for an arbitrary number of phases, again without area constraint, and refor- 
mulate the algorithm in a vector-type setting. In the third step we finally 
design the method for multiphase area-preserving motion. 

2.2.1. Two-phase motion without area constraint 

We describe the BMO algorithm for the case when only two phases are 
present. This algorithm works in any space dimension. Given an initial 
interface 7, we take P to be the region enclosed by this interface (and possibly 
by the boundary of the region Vt where the motion is considered) and define 
its characteristic function x 



of the solution to the heat equation at time At with initial datum x- The 
two-phase algorithm thus reads as follows: 
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1. Given a region P, set x to be its characteristic function. 

2. Solve the heat equation with initial condition x'- 

ut{t,x) = Au{t,x) for {t,x) e (0, At] x Q, 
Ou 

— = on (0, Atl X on, 

on 

u{0,x) = x{^) inQ. 



Update X as the |-level set of u{At, 



1 if u{At,x) > |, 
if ulAt,x) < i. 

The evolved interface is now the boundary of the set {x G Q; x{^) = !}■ 
4. Go back to step [2] to proceed with the computation for the next time 
step. 

It has been rigorously shown, in a general framework including topological 
ch ang es, that this algorithm converges to motion by mean curvature as At — >■ 

(iirsQ- 

Here we remark that the Neumann boundary condition in the diffusion 
step guarantees that the interface will touch the boundary of Q with right 
angle. Other boundary conditions, such as Dirichlet conditions pinning the 
interface at the boundary, may be used according to necessity. 



2.2.2. Multiphase motion without area constraint 

We next address the case of multiple phases. The idea of sharpening sep- 
arately diffused characteristic functions (one for each phase) was introduced 
in j22|. The algorithm is as follows. 

1. Given regions Pi,i = 1, . . . ,k, set Xi to be the characteristic function 
of P.. 

2. For i = l,...,k, obtain Ui{At,x) by solving the heat equation with 
initial condition Xi up to time At. 

3. Update Xj ^ the characteristic function of the set where Uj has the 
largest value amongst all solutions 




if Uj{At,x) > Ui{At,x) Wi 7^ j, 
otherwise. 



The new interfaces are the boundaries of the sets {x G f2; Xi{x) = !}• 
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4. Go back to step H] to proceed with the computation for the next time 
step. 



The above algorithm can be reformulated to obtain an equivalent algo- 
rithm using a single vector-valued heat equation. This is essential for imple- 
menting the area constraint and for dealing with more general motions. 

We prepare k reference unit vectors Pi,i = 1, . . . ,k, of dimension k — 1, 
each corresponding to a phase Pi. They are defined as the vectors pointing 
from the centroid of a standard fc-simplex to its vertices (cf. figure A. 12 in 
the appendix). Hence, there are infinitely many possible A;-tuples but the 
relative distributions of the vectors are identical. See Appendix A for a 
simple way to construct these vectors. 

Using the vectors p^, the multiphase algorithm can be written as follows: 

1. Given regions Pi,i = 1, . . . , k, set 

Uq{x) = Pj for X G Pi. 

2. Solve the vector- valued heat equation with initial condition uq: 

utit.x) = Au{t,x) for (t, x) G (0, At] X (3) 
du 

— (t,a;) = on(0,At]x on, 

on 

u{0,x) = uo{x) inQ. 

3. Update no by identifying the reference vector which is closest to the 
solution u{At, x): 

Uo{x) = Pj, where Pj ■ u{At,x) = rnax^Pj ■ u{At,x). (4) 

This redistribution of reference vectors determines the approximate new 
phase regions after time At. 

4. Go back to step [2] to proceed with the computation for the next time 
step. 

The equivalence of our algorithm with the original BMO can be shown 
by considering the functions 

Wi{t, x) = ^-j-^ (u{t, x) ■Pi + -r-^] , i = l,...,k. 
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Since u is the solution to the hnear vector- valued heat equation, we have 
{wi)t = Awi in (0,T] X n, 

^ = on (o,r] X dn. 

on 

Moreover, since from the construction of Pj it holds that ■ Pj — 1/(1 — k) 
for i j, we can readily check the identity 



where Xi is the characteristic function of i-th phase region. It follows that 
Wi is identical to the solution Ui of the scalar heat equation from the original 
algorithm for each i — 1, . . . ,k. Prom the definition of Wi it is immediate 
that 

u{x,t) • Pj > u{x,t) ■ Pj <^ Wi{x,t) > Wj{x,t) <^ Ui{x,t) > Uj{x,t), 

which proves the equivalence of both algorithms. 

Remark. The reference vectors are related to the position of wells in the 
phase field approach. Indeed, the idea of the BMO algorithm for the case 
of two phases originated in a simple splitting scheme for the singularly per- 
turbed reaction-diffusion equation 



where 1^ is a double-well potential. Here, the sphtting scheme consists of two 
steps. The first step solves the heat equation Ut = eAu, which corresponds 
to the diffusion step of the BMO algorithm, and the second solves Ut = 
— ^W'{u). This corresponds to the thresholding step if the equation is solved 
for sufficiently long time. Here we can see that the thresholding values and 
1 in the BMO algorithm correspond to the positions of the two wells of W. 
Accordingly, if we want to calculate three-phase motion, we can look at the 
problem from the viewpoint of constructing a suitable well-type potential. 
A potential with three wells at different positions along a scalar axis would 
obviously yield incorrect results, because the strength of the wells would not 
be equivalent. Therefore, we have to increase the number of variables for the 
potential and construct the wells in a symmetric way. The reference vectors 
introduced above then give the coordinates of the positions of the wells. 




ut^eAu- lW'{u), 
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2.2.3. Multiphase motion with area constraint 



The paper [23| presents an area-preserving BMO method for 2 phases. 
The authors adjust the height at which the thresholding occurs in such a 
way that the resulting area of the level set is preserved. Such a level set is 
guaranteed to exist by the maximum principle and one can compute that the 
thresholding height has to be changed from | to 




where Ka is the average mean curvature along the interface. 

However, when three or more phases are present, each interface has a 
different average mean curvature and the thresholding heights become dif- 
ferent. One could try to use a different thresholding height for each phase 
but then the global interaction would be ignored and the phases may overlap 
or create vacuums, especially when they are initially touching. Therefore, 
we suggest a different approach incorporating the area constraint into the 
diffusion process. In this method, a heat equation with a nonlinear source 
term expressing the area preservation of level sets is solved. Subsequently, 
the solution is sharpened at a fixed height. 

The mentioned nonlinear heat equation corresponds to the gradient flow 
of the functional 

J(u) = [ \Vu\'^dx 

in the constrained set 



n 



K, = ^^u E H^{^l;R^ / X{uixyp,>uixyp^ vj} dx = Ai for i = 1, . . . , A; - l|, 

where Ai denotes the given area of phase Pi. For simplicity, we consider only 
the case when one type of interface divides one type of phase into multiple 
regions, like in a soap froth. When multiple phases and interfaces with 
nonuniform properties are present, a different form of the energy J has to be 



adopted ^ 



Let us consider the two-dimensional two-phase case as in figure 1 to un- 
derstand the meaning of this gradient flow. Our argument is formal and 
we simplify the exposition by discretizing time and adopting the numerical 



approach from the next section, which is based on the method of Rothe [29 
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Let the curve 7 denote the |-level set of a scalar function u and call the 
region enclosed by this curve P. We assume that an initial function u^, is 
given and we search for the minimizer of the integral 

J(u) = [ f + \VuA dx (5) 



h 



under the constraints 



m(7(s)) = ^, se[a,6], (6) 
meas(P) = A. (7) 

Here h denotes the length of the discrete time step and A is the required 
area. 

Note that by perturbing function u in regions away from the ^-level set, 
one can readily deduce that the minimizer u satisfies, in a weak sense , the 
following 

^i—^-Au = inPU(^]\P), (8) 

^ = ondn. 
on 

Now consider a perturbation of u of the form u + Su, which is allowed 
to affect the |-level set. The corresponding change in the level curve can 
be written in the form 0(7 + ^7), ^7 = (571,572), where a is a constant 
depending on ^7 and whose role is to adjust the area to the correct value. 
Because of constraints and ([7]), we cannot choose 6u arbitrarily. However, 
we can select an arbitrary ^7, find an appropriate constant a, and use the 
corresponding 6u. 

We compute the variation of functional J remembering that u may not 
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be smooth across the interface 7: 

J{u + 5u) — J{u) 
' u — u 



^ -5u + VuV{5u) ] dx 



h 



5u + VuV5u ] dx + [ — ; — -6u + VuVSu ) dx 
J Jn\p \ h J 



Am ) 5udx+ I ^^^^5updl 



h / dn 



u — 



Au\5udx— / — —^5uQ\pdl+ I —Sudl. 
n\p \ h J dn Jq^ dn 

Here wp denotes the value of a function w taken as a hmit from inside of P. 
The symbol if n\p has the analogous meaning. 
Due to dH]), on the interface we have 

dup 9uQ\p \ 

— — 6up — SuQ\p ] dl = for all admissible 6u. 

on on J 

We next to rewrite this identity in terms of the arbitrary perturbation 6'y. 
To this end, we use the fact that the phase area is preserved to obtain 

^c? J [(71 + 571) (72 + '^72) - (72 + S-f2){l'i + h'l)] ds = A, 
which yields 



S'j ■ ndl ^ — A. 

a"' 

Here, the symbol ~ means that the equation holds up to second order in 
terms of ^7 or Su. From this we can compute the value of a — 1 which will 
be needed later: ^ 

a — 1 ~ / ^7 ■ n dl. (9) 

Condition ([6]) means 

(m + (5m)(q;(7 + 57)) — ^(7) = 0. 



16 



Using a Taylor expansion and we obtain that, on 7, one has 
~ Vwp ■ ((a — 1)7 + 57) + 5up 

~ Vup ■ ^57 — ^ J 6j ■ ndl^ + 6up 



(10) 



n\p 



(57- 



2A 



5'y ■ ndl) + 5u. 



''Q\P- 



Expressing 6u from (JTOll as 



5m 



9u/ 
(9m 



n 



n\p 



n\p 



dn 



n 



57 



2A 



6^ ■ ndl 



we obtain 



,duc 



'7 L 

Noting that 



dn 



n 



n 



5'j — 



1_ 
2A 



5-i-ndl d/ = V57. 



7 



^7 — — / 5'^ ■ ndl] dl = 0, 



we arrive at the interface condition 



dup 
dn 



du 



Q\P 



dn 



A = const 



on 7. 



Similar calculation can be carried out in the multiphase setting, see Appendix C 
for more details. 

In view of the derived condition on 7 and the results regarding two-phase 
free boundary problems in [30], we can reformulate the variational problem 
in the following way: Find a minimizer u\ G (Q) of the functional 



\u — u 



h 



h iVwp + Ax„>i ) dx. 



(11) 



According to the results in [30|, we can expect that this minimization problem 
is in a sense equivalent to looking for a weak solution of the problem 



Ut — Au = in (0,T] x Q, 

on dQ, 



[12) 



du ^ 
dn 
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where /x is a Radon measure given by 

x) = \{t)V} LaP(t), P{t) = {x; u{t, x) > i}, 

for a suitable space-independent function A. Here the symbol "H^ means the 
one-dimensional Hausdorff measure. This type of problem is known as the 



two-phase parabolic free boundary problem [31|, l32|, |33 . 

Having deduced the partial differential equation corresponding to the 
constrained gradient flow, let us return to the considerations concerning the 
BMO method. Formal calculation shows that the BMO method using the 
parabolic equation f ll2p instead of the heat equation evolves interfaces with 
normal velocity equal to minus mean curvature, plus a space-independent 
term: 

w = -K + 2A(0) + 0(t), t-^0 + . 

Thus, if A is chosen appropriately, the area-preserving mean curvature flow 
is realized. Since the derivation of the above relation is rather technical, we 



present it in Appendix B 



Minimization of fill I) is much easier to implement than the constrained 
minimization or the partial differential equation (11 2p . In each time step of 
the BMO algorithm, if for some A we obtain that the area of {x; ux{x) > |} 
is equal to the given value A, then ux becomes a solution to our original 
problem. However, it is not clear how to calculate such a A or even if such 
a function exists. For the two-phase case we see that A should be half of 
the average mean curvature, but determining the Lagrange multiplier for the 
multiphase case is complicated. Moreover, one of the advantages of the BMO 
approach is that it does not require the computation of curvature values, so 
we want to avoid the direct calculation of A. 



Therefore, based on the results in [Sj] we approach ([5]) by using the 



method of penalization and consider minimization of the functional 

Je{u)= I ( '"^ + |Vn|^) dx + P,{\{u > III - A), 

where e is a small positive number and the penalty function is defined by 

Pe{s) = I '''' ' - (13) 

^ ' \ -es if s > 0. ^ ^ 

Penalization techniques of this type for stationary problems were analyzed 
in SJ, |35|, |36| and other works. The general conclusion of these treatments is 
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that the penahzed solutions converge as £ —t- to a solution of the correspond- 
ing PDE or of the original constrained minimization problem. Moreover, it 
was discovered that, in order to obtain the solution of the original problem, 
it is not necessary to carry out the limit e — )■ 0, since the exact minimizer 
is obtained for some sufficiently small but positive value of e. Although the 
analysis of the time-dependent problem has not yet been addressed, this ob- 
servation is essential to justify the robustness of the penalty method, at least 
in the stationary case. 



Remark. A hint at another approach to approximating the constrained min- 
imization can be obtained from the works developing the theory of singular 



perturbations [3l|, |32|, |33|]. The solution is realized as a limit of solutions to 
heat equations singularly perturbed by a nonlinear source term. The results 
also cover the time-dependent problem but are at the moment limited to free 
boundary conditions independent of the solution. The application of this 
approach to problems of the type f|T2l) . where A is a nonlocal term depending 
on the solution (such as the average mean curvature of level sets), is an open 
problem, which we would like to address in the future. 



In closing this section we note that, except for the paper [36|], the con- 
strained multiphase (vector-type) problem is almost unexplored. We provide 
here a formal analysis of the convergence of the BMO algorithm for the con- 
strained multiphase problem. Although this analysis is an important part of 



the present work, we defer it to [Appendix C| in order to make the main text 
less tortuous. 



3. Numerical computation 

In this section, we present a discrete algorithm for realizing constrained 
multiphase mean curvature flow and give the details of its numerical imple- 
mentation. 

3.1. The numerical algorithm 

In treating the multiphase motions, the reformulated BMO process, stated 
in terms of a vector-valued heat equation as in section \2.2.2\ is approxi- 
mated by use of a minimizing movement. A discretization in time is used 
to build approximate solutions by successively minimizing time-independent 
functionals; hence this setting conveniently allows one to include constraints 
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via penalization. Each minimizer then corresponds to the solution of a 
vector-valued elliptic problem with Lagrange multipliers appearing on the 
free boundaries. 

Namely, the heat equation of the BMO algorithm is solved by means of a 



vector- type discrete Morse flow (DMF) [29|,l37j. In the unconstrained case, at 
each step of the BMO algorithm we have to solve equation for a time At. 
Given regions Pi, i = 1, . . . , k, we start the DMF by constructing a function 
Wo such that Uo{x) = if x E Pi, and set Wq = Uq. We choose a large 
positive integer K which determines the time step of the flow, h = At/K. 
Then, to obtain the approximate interface position after time At, for n = 
1, ...,K we inductively minimize the following functional over H^{Q; R'^^^): 

^•■'"'> = X( 2ft + VJ''" 

To evolve the interface for a time T, our method takes M = T/At and 
repeats the following for m = 1, M: 

1. Set Wq = Wm-l- 

2. For n = 1,...,K, compute tu„ to be the minimizer of J'n{w) over 

3. Obtain Um by thresholding: 

Umix) = Pj, where Pj ■ wk{x) = max^p^ ■ wk{x). 

The sequence of functions {um}m=o then gives an approximation to the (un- 
constrained) multiphase motion. Figure 3 shows the basic characteristics of 
the method for a four-phase problem, but see section H] for a clarification of 
the initial condition. 

In treating area-constrained motions, one can see that the minimization 
aspect of our reformulated algorithm should allow the inclusion of area con- 
straints via penalization. For example, denoting the prescribed area of region 
Pi by Ai, the energy functional can be modified: 

k 

J'niw) = Mw) + -Y1 1^' - "^eas(i^r)l'- (15) 

i=l 

Here e > is a small penalty parameter and the areas corresponding to w 
are obtained from the sets 

PI" = {xe n; w{x) ■ Pi > w{x) ■ p. Vj}. 
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Here we note that the penalty fllSp . which is used in our numerical com- 
putations, is slightly different from the theoretical form f ll3p . However, we 
prefer this form, since it is more simple and gives satisfactory results. 




Figure 3: (Left) The initial condition. (Center) An instant in time of the solution to the 
vector-valued heat equation. (Right) The corresponding interfaces. 



3.2. Implementation of the method 

The numerical implementation of our method uses the finite element 
method to approximate the functional values (fTSll . and minimizers are found 
by gradient descent. The domain is triangulated into a finite number of el- 
ements, over which we assume that the function is continuous and a linear 
interpolation between node vectors. The solution to the vector-type equa- 
tion (IC.8|) is thus approximated by successive minimizations of f lTSj) until 
arriving at the thresholding step. As noted in the introduction, and which is 
well-documented, simple thresholding by (jl]) is known to inhibit the motions 
obtained by computing with the BMO algorithm (see 381] ). We now briefly 
explain the troubles which may occur. 

Before thresholding, the interface is a level set of the finite element solu- 
tion to the heat equation, and, in the volume-preserving case, approximately 
satisfies the prescribed area constraint. However, applying the original for- 
mulation of the thresholding (|1]) at the nodes of the mesh would then alter 
the position of the interface. This causes two difficulties, the most signifi- 
cant being that, upon proceeding to the next minimizer, the interface may 
fail to move (thus becoming stationary). That is, since each evolution is ob- 
tained via a heat equation, the diffusion process must proceed long enough 
so that the grid resolution resolves the movement of each interface across 
the elements. On the other hand, if the diffusion proceeds too long, the 
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approximation of the mean curvature flow looses its accuracy; hence certain 
configurations do not permit a suitable time step. Additionally, due to the 
constraints on the area of each phase, the normal velocity of the interfaces 
tends to be much slower than that of the unconstrained motion, especially 
near the stable state. Therefore this issue is particularly relevant to our cur- 
rent problem. The second issue is that the enclosed areas cannot be preserved 
with sufficient precision after this thresholding. 

We are able to overcome these issues by use of the following. Just prior 
to the thresholding step (i.e., after obtaining wk), we indicate the elements 
that span interfaces by e* , record the interfacial geometry, and then threshold 
by dl]). Upon the next minimization, which is the first minimization of the 
next BMO step, whenever we come to an indicated element, we recall the 
geometry of the interfaces and compute the value of the functional f|T5|) by 
means of a triangulation of the element. In particular, the value over an 
indicated element is obtained by the values over a set of convex polygons, 
each denoted by Ri. These regions are determined by the element nodes and 
the intersection of the recorded interfaces with the element edges (see figure 
4): 

Ri={xee*; Wk{x) ■ Pi> Wk{x) ■ Pi VI ^i}. 

For a candidate minimizer u, the value of the discretized term in the func- 
tional dH]) is then computed: 

and the contributed for the penalty terms 

meas{e* fl Pi) = meas{Ri). 

are accumulated. By such an approach, we are able to realize interfacial 
motions whose configuration and precision of enclosed area are not altered by 
truncations. Moreover, this approach allows one to alleviate the restrictions 
on the time and space discretization of the standard BMO algorithm (see the 
analysis in section I4.1.ip . 

4. Numerical tests 

This section presents numerical analysis and a number of numerical ex- 
amples of the application of our method to area-constrained flows. We use Pi 
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Figure 4: A stable solution, a junction close-up, the partition scheme for the vectors 
(figures use the real data). 



Lagrange finite elements in each of tlie computations. Tlius, after assigning 
tfie appropriate vector G R'^"^ to every node of tfie mesh (see figure A. 12 
for low-dimensional visualizations of these vectors), we note the jagged shape 
of each initial condition. Although curvatures are not defined for such initial 
conditions, our diffusion-based algorithm is able to handle their evolution 
without trouble. 

The physical interpretations follows. We configure a given number 

of bubbles into the shapes shown in the figures and then let them evolve. 
As we assume that there are no inertial forces, the evolution by the steepest 
descent (with area preservation) of the length energy ([2]) can thus be thought 
of as expressing the slow movement of the bubbles. By examining the data, 
we note that the so-called symmetric Herring condition (junctions meeting at 
120°) appears to be satisfied at junctions (see figure 4 for a close-up inspection 
of one such junction). 

We again mention that the process described in ( [T6l) is essential in com- 
puting these motions. Indeed, as the area-preserving interfacial evolution 
tends to be slower than motion by mean curvature flow, the well-known time 



and grid spacing restrictions of the BMO become particularly relevant [39 
Nevertheless, recalling the interfacial geometry after thresholding allows us 
to avoid such complications, and can also be used for the non-constrained 
BMO with the same result. That is, our method also works for large ratios 
of grid and time step sizes, for which the original BMO becomes stationary. 
Moreover, formal analysis and numerical tests show that this approach does 
not alter the characteristics of the target motion. 
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It is also worth stating that another technique for handhng the restriction 
of the BMO algorithm on non-adaptive meshes is through the use of signed 



distance functions. This method was developed in |38[, where it was shown 
that it gives satisfactory results and we would like to extend the method used 
there to the multiphase constrained case. 

4.I. Convergence analysis 
4-1.1. Analysis of error 

This section examines the behavior of our method in comparison to the 
standard BMO. We refer to the standard method as BMO, and to our own 
algorithm (which utilizes the thresholding from section as BMO*. 

We examine the error of our method when applied to a simple test prob- 
lem. By symmetry, a circle of initial radius ro which is evolving by mean 
curvature flow remains a circle whose radius r{t) satisfies the following ordi- 
nary differential equation: 

r' = -- (solution: rit) = (rg - 2ty/^). 
r 

With an initial condition obtained from the target radius ro = 0.35, we 
vary grid and BMO time step sizes, and compute the solution by use of the 
BMO and BMO* algorithms. The output of the program is a list of interface 
nodes {{P^ , P^)}i. Circles were fitted to these points at each time level by 
minimization of the functional 



(^{p,' - c^y + {P! - cyf - r^y 



with respect to the centre coordinates {C^, C^) and radii r. We measure the 
error of the method by the time-average of the absolute difference between 
the radius of the fitted circle and the exact radius: 



1 ^ 

7$^|rfit(tz)-r(tz)|. 



Here L denotes the number of time steps until the radius is zero. The error 
table for the standard BMO algorithm is as follows: 
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res(space\time) 


2 


4 


8 


16 


32 


64 


128 


256 


5x5 


0.0483 


0.0231 














10 X 10 


0.0572 


0.0176 














20 X 20 


0.0047 


0.0018 


0.0166 


0.0070 










40 X 40 


0.0044 


0.0039 


0.0034 


0.0036 










80 X 80 


0.0056 


0.0032 


0.0027 


0.0101 


0.0042 


0.0223 






160 X 160 


0.0060 


0.0031 


0.0028 


0.0073 


0.0107 


0.0039 


0.0040 


0.0397 



We confirm that when the time step decreases to a certain size (relative 
to the gird), the BMO method halts (this is indicated by the symbol "-"). 
The critical ratio of space to time step size is approximately 20. 

On the other hand, by computing the same evolutions with the BMO* 
algorithm, we find that it is able to deal with a much wider range of parame- 
ters (critical ratio around 400). The initial condition is an approximation to 
a circle of radius ro = 0.35 obtained from the area-preserving BMO* method 
with penalty parameter e = 10~^. The error table is as follows: 



res(space\time) 


2 


4 


8 


16 


32 


64 


128 


256 


5x5 


0.0276 


0.0228 


0.0258 


0.0381 


0.0440 


0.0368 






10 X 10 


0.0121 


0.0146 


0.0076 


0.0165 


0.0201 


0.0103 


0.1210 




20 X 20 


0.0044 


0.0038 


0.0033 


0.0043 


0.0081 


0.0139 


0.0080 


0.0694 


40 X 40 


0.0045 


0.0035 


0.0024 


0.0024 


0.0028 


0.0007 


0.0056 


0.0096 


80 X 80 


0.0053 


0.0032 


0.0025 


0.0042 


0.0070 


0.0073 


0.0056 


0.0016 


160 X 160 


0.0059 


0.0031 


0.0026 


0.0049 


0.0085 


0.0097 


0.0097 


0.0088 



Overall, we see that both algorithms approximate the solution, and that 
the additional partitioning of the BMO* algorithm is beneficial. For finer 
meshes, both algorithms show a tendency of stagnating errors when the space 
mesh is refined. This fact could be attributed to the properties of the DMF 
scheme. However, we note that the point of this partitioning is not to improve 
the error per se, but to relax the time and grid restrictions inherent to the 
standard BMO. 

4.1.2. Analysis of penalty parameter 

Here we perform an error analysis for the two-phase area-preserving case. 
Since the velocity of the area-preserving motions tends to be much slower 
than in the non-constrained case, the use of the BMO* should be preferred 
over the standard BMO. 

Wc take two non- intersecting circles of radii ra — 0.1996 and ri, — 0.1384 
and identify them as the same phase. Then the area preservation condition 
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implies that the radius of the larger circle will grow as the smaller circle's 
radius shrinks. The evolution of the radii follows the equations 

-1 2 

— + 

ri ri + r2 

-1 2 

— + . 

We use a numerical method to compute a precise approximation to the so- 
lution of the above system and compare the result to that obtained by our 
penalty method. In the BMO* method, the domain [0, 1] x [0, 1] is triangu- 
lated into 14536 elements, At = 2.5 x 10""^, and K = 30. The results are 
ploted in figure 5. 

Weak penalties with e = 10°, 10-\ lO'^ give almost the same result, 
and the larger circle shrinks slightly although it should grow. For increas- 
ing strength of the penalty the results approach the correct solution with 
the larger circle growing as the smaller shrinks. Finally, penalties e = 
10^^, 10"^, 10~^ give almost identical results, which are only slightly differ- 
ent from the exact solution. This finding agrees with the theoretical predic- 
tion mentioned at the end of section I2.2.3j namely that we should obtain 
the exact solution for sufficiently large penalties. The solutions for large 
penalties indeed do not change, and their slight deviation from the exact 
solution is caused by the discretization. In conclusion we can say that the 
penalty method behaves well, since for each grid resolution there exists a 
range of penalty coefficients for which the solution is independent of the 
penalty strength and appropriately approximates the exact solution. 

4.1.3. Analysis of the multiphase area-preserving algorithm 

In order to test the performance of our multiphase area-preserving algo- 
rithm, we compute the stationary solution corresponding to two 2-dimensio- 
nal soap bubbles attached to a wall. It can be rigorously shown that the 
steady shape is composed of three circular arcs meeting with 120° angles at 
the triple junction and with 90° angles at the walls. Moreover, the radii of 
the arcs satisfy the well-known condition 

1 1 _ 1 

ri r2 ri2 ' 
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time 



Figure 5: Evolution of the radii for penalty parameters e varying from 10° to 10 ^ (black 
curves, ordered from left to right). The red curve shows the exact solution. 



where ri and r2 are the radii of the bubbles and is the radius of their 
common arc. When the initial volumes of the bubbles are given, it is possible 
to compute the above radii analytically. 

This analytical solution is compared to the numerical solution in figure 6. 
The numerical solution was obtained by running the area-preserving three- 
phase method for sufficiently long time, until the interfaces stopped moving. 
The domain [0, 1] x [0, 1] was triangulated into 14536 elements. At = 2.5 x 
10~^, K = 30, e = 10~^, and the fitting method from section 14.1.11 was 
used to obtain the circle radii. Although the area differs slightly from the 
prescribed value, the configuration of the numerical solution, including the 
triple junction, agrees well with the analytic solution relative to the resolution 
of the grid. 

4-2. Multiphase flow 

Here we present two examples of multiphase mean curvature fiows. 

4-. 2.1. A triple bubble 

We begin by examining the motion of four phases, where the three bubbles 
initially touch each other (figure 7). The area of each phase is different, and 
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0.1 0.3 0.5 0.7 0.9 0.54 0.56 0.58 0.6 0.62 

Figure 6: The stationary configuration for two bubbles attached to a waU: (Left) Exact 
solution (blue) and numerical solution (red) - circles obtained by least squares fitting are 
plotted. (Right) A close-up view of the triple junction showing the exact solution (blue), 
fitted numerical solution (red) and original numerical solution (black with dots). 

a corresponding triple bubble is obtained as the stable configuration for large 
times. 

Here the domain Q = [0, 1] x [0, 1] is triangulated into approximately 5600 
elements, At = 5 x 10~^, and K = 10 . The mesh is nearly uniform so that 
most elements have area approximately equal to 1.8 x lO"''. A penalty of the 
form shown in (fT5|) is added for each phase and its parameter is e = 10~^. 
The prescribed volumes were maintained to within an absolute error of 10~^, 
even during the dynamic portion of the movement. 




Figure 7: (Left) The initial condition. (Center) Evolution after a short time. (Right) The 
stable configuration. 
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4-2.2. A nine phase flow 

In the present computation, nine phases are positioned throughout the 
domain. After assigning the appropriate vector to each node location and 
detecting the interfaces, we have the initial interfaces as shown in figure 8. 
The bubbles cling to the lower boundary and, eventually, the top most bubble 
slides itself in between two others. This shows that the method can naturally 
handle topological changes. 

The domain Q = [0, 1] x [0, 1] is triangulated into approximately 9600 
elements. At = 3 x 10~^, and K = 30. The mesh is nearly uniform so that 
most elements have area approximately equal to 10~^. A penalty of the form 
shown in ( 1TB]) is again added for each phase and its parameter is e = 10~^. 
The absolute errors in the areas were able to be maintained within 10~^. 




Figure 8: The initial condition, evolution after SOOAt, and the stable configuration. 



4.3. Interaction of interfaces 

The proposed method can also handle interactions of interfaces, such as 
merging and attaching. 

4.3.1. Two-phase coalescence 

We first consider a computation involving two phases. One phase is ini- 
tially separated into two distinct parts and is configured in such a way that, 
in the course of the interface evolution, the two parts eventually come into 
contact with each other. The parameters are as in section 14.2.11 

When the two parts touch, a topological change occurs, and the evolution 
continues until reaching a stable configuration (a circle); see figure 9. We note 
that, due to the diffusion process of our method, the initially separate phases 
attract each other slightly. Nevertheless, this attraction decreases extremely 



29 



fast with the distance of the phases. However, for this reason, capturing 
the precise behavior at the moment of the topological change is difficult. Of 
course, this unwanted attraction can be reduced by refining the time step 
and grid resolution. The prescribed volumes were maintained to within an 
absolute error of 10~^. 







0.25 0.5 0.75 1 



Figure 9: Merging of two regions corresponding to the same phase. 



4-3.2. Three phase coalescence 

In this computation, we place three phases throughout the domain. Two 
are configured to be initially separate, but so that they eventually touch. As 
the areas of the phases are different, the final steady state solution shows 
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the shape of a non-symmetric double-bubble; see figure 10. The setting of 
parameters is identical to that of section 14.2. 1[ Absolute errors in the areas 
are approximately 10~^. 




0.25 0.5 0.75 





0.25 0.6 0.76 1 




025 0.5 0.75 



Figure 10: Attaching of two regions corresponding to different phases. 



4.4- additional transport term 

Here we deal with a generalization of the constrained mean curvature flow 
to include a transport term. For simplicity, we shall focus on the two-phase 
case which can be described by a scalar model. 

Let us consider the partial differential equation (compare to (IT^ ). 



Ut = Au 



f 



'Ant 



(17) 



where / is a specified function of space. It can be shown in a fashion similar 



to Appendix B that the application of the BMO process to (fT7|) leads to the 
motion of interfaces with normal velocity 



V = — K — / + A, 
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where k is the mean curvature and A is a function of time that guarantees 
the preservation of area. 

In the numerical computations, one uses the method described in section 
13.11 minimizing the following penalized functional over H^{VL): 

where A is the enclosed area of the bubble which should be preserved over 
time and A{w) is the measure of the set {w > 1/2}. 

We apply the explained method to carry out a simple two-dimensional 
simulation of gas bubbles rising from the bottom of a container filled with 
a viscous hquid (see also 4o|). In this case we set / = /Sy, where y is the 



coordinate direction of gravity and /3 is a constant expressing buoyancy. We 
consider the case of a bubble having the shape of a partial ellipse and initially 
positioned at the bottom of the container. Figure 11 shows the evolution at 
four distinct times for two different initial shapes. The results were obtained 
using the parameters At = 10-^ K = 20, e = 0.001, and (3 = 20.5. 

We compute under Neumann boundary conditions which means that the 
bubble will always touch the bottom with right angle. The motion is a result 
of the balance between the buoyant force pushing the bubble upwards and 
the surface tension force pressing the bubble towards the bottom. For the 
bubble on the left, the adhesive force is prevalent so it attaches itself to the 
boundary and then comes to a rest in a stationary shape. The bubble on the 
right has a shape for which the buoyant force wins and the bubble detaches 
itself from the boundary. After detaching, the bubble becomes circular and 
continues moving upward. 




Figure 11: Volume-preserving mean curvature flow with buoyancy - comparison of evolu- 
tion for different initial shapes. 
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The possibility of a straightforward inclusion of a transport term into the 
numerical algorithm as explained above is extremely important in applica- 
tions where two or more phases interact at the interface, such as through 
buoyancy in our simple example. We believe that it is possible to proceed 
in this direction and design effective algorithms for various interaction prob- 
lems. 



5. Conclusion 

We developed a method for approximating constrained multiphase cur- 
vature-driven motions. Our method is based on the BMO algorithm, which 
was reformulated in terms of a vector-valued heat equation. We derived the 
nonlinear PDE which governs the corresponding constrained evolutions and 
used it to formally prove convergence in the multiphase setting. 

The vector-valued BMO algorithm was implemented employing the dis- 
crete Morse flow and it was found that the variational nature of this approach 
allows one to consider constraints via additional penalty terms. Using this 
idea we were able to realize multiphase area-preserving interfacial motions 
that are free of the defects of other methods. 

By detecting the precise locations of the interfaces we are able to com- 
pute the area of each phase at a high precision and thus impose the area 
constraints. This geometric information is also retained after the thresh- 
olding step, which was found to alleviate the BMO's standard restrictions 
on the spatial and time-step resolutions, for both the constrained and non- 
constrained problems. 

In closing, we remark that it would be interesting to investigate our algo- 
rithm in relation to the recent threshold dynamics utilizing signed distances 



functions [38|, and to consider its position in apphcations. 



Appendix A. Construction of reference vectors 

For the sake of completeness, we include a method for constructing the 
vectors corresponding to the regular simplexes. 

The computation of the reference vectors can be done by first considering 
the standard (A; — l)-simplex in and rotating its vertex vectors in a suitable 
way. Particularly, the standard {k — l)-simplex is given by 

k 

Sk-i = {{xi, X2, ■ ■ ■ , Xk) G M^; Xi = 1, and > for i = 1, . . . , k}. 

i=l 
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Its vertices have the coordinates 

(1,0,. ..,0), (0,1,0,... ,0),..., (0,...,0,1) 

and if we translate the simplex in such a way that its centroid lies in the 
origin, the vertices will be 

pl = i(_i,...,_i,/,_i). 

We fix an orthonormal basis for the (/c— l)-dimensional hyperplane containing 
the simplex. A convenient way is to take the first translated vertex above 
as the first basis vector (after normalization) and construct the remaining 
vectors as follows: 

q'l = , ^ (fc- 1,-1,-1,-1 -1) 

«' = TlriKS^C.''.*- 3.-1. ■■■.-!) 

9^1 = ;^(0,..., 0,1,-1). 

Denoting the matrix having q'l, . . . , ql_i as its rows, we obtain the ref- 
erence vectors as the normalized projection of = l,...,k, into this 
orthonormal system, i.e. 

Appendix B. Formal proof of convergence of the modified BMO 
algorithm 

We show that the |-level set of the solution to the problem 

ut = Au + n in(0,T)xM2, 
u{0,x) = Xp{o){x) inM^ 
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Figure A. 12: (Left) The 2-phase regular simplex. (Center) The 3-phase regular simplex. 
(Right) The 4-phase regular simplex. 

evolves according to its curvature k, plus a constant factor A, plus an error 
term which approaches as T — >■ 0. Here P(0) is a given initial region, P{t) 
denotes {x : u{x, t) > |}, and is a Radon measure given by 

for a suitable function A. We assume that A(0) can be defined so that A is a 
smooth function in [0,T]. 
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Figure B.13: Configuration of the interface in the proof of BMO convergence. 

Wc consider the coordinate system as in figure B.13, where the point 
(0, 0) lies on dP{0) and the outer normal n to P{0) at this point is the vector 
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(0, 1). We assume that inside the cube Q = {{xi,X2) : |a;i| < 1, \x2\ < 1} the 
boundary of -P(O) is given by the graph of a function 7 : (—1, 1) — )■ (—1, 1), 
satisfying 

|7'(2;i)|<l, XiG(-l,l), 7(0)=7'(0) = 0. 

Then 7"(0) is equal to minus the curvature k of dP{0) at the point (0,0). 
Moreover, we assume that for a sufficiently short period of time, dP{t) can 
be written as the graph of a function 7(t, xi) in the same coordinate system. 
Let V be the normal velocity of P(0), i.e., 

vtn e dP{t) or u{t, 0, vt) = i. 

Writing the explicit solution to the heat equation with a source term, we 
have that u{t, 0, vt) is equal to 

1 f 4+(=^2-vtf ft r 1 ^j+(^2-^o^ 

e 4* ax + — -e ^(f^^) fi{s,xj ax as. 

4vrt Jp(j) Jo yiR2 47r(t - s) 



By the results of |26l. |41[|. the first integral on the right-hand side is equal to 

Let us denote the second integral on the right-hand side by /. Our goal is to 
show that 

/ = ^A(0) + O(t3/2), t-^0 + . 



We split the integration domain into two parts as follows: 

/■* r x(s) 4+(-2-tf 

Jo JdPis)Mt-s) 

f X(s) 4 + (x2~vt)^ 

/ ^ ^ e dS{x)ds 

JdP{s)nQ 47r(t - s) 



/ ' , e ^(^^^dSix) ds 
JdPis)\Q 47r(t - s) 



h + h 
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First we show that the second integral I2 is exponentially small: 

r f \m\ ^ ■ 

'0 JdPis)\Q 4vr(t - s) 



I/2I < 



-eW=^dS{x) ds 



< max \dP(s) 
selo,t] 



Jo 



* As 
47r(t - s) 



C \\{t - i^)\-e-'^ da 

Jl/{4t) ^ 



< 4Ctmax|A(s)| / 



e-'^da 



= Cte-^t 
= C>(e-"/*), 

for some a > 0. In the above calculations we use the change of variable 
a — l/A{t — s) and we assume that dP{s) has finite length for s e [0, T] and 
that A(s) is bounded for s e [0,T]. 

In the case of Ii we use the fact that the boundary of the domain is a 
graph and employ Green's theorem. We use the notation p = y^t^^ in 
order to simplify the formulae. Further, we perform a change of variables 
zi — xi/ p, Z2 = {x2 — vt)/p to obtain the value of Ii as 



'0 J-iJ-00 



dxn 



Arr{t - s 



x^-\-(x2—t^t)^ 

-e 



^1 + 7'(s, a;i)2 dx2 dxi ds 



^ jHe ^<*-> ' \^l + i{s,xiYdx2dxids 



t l-l/p 



7{s,p2]^ ) — vt 



A(s)„_ 



4+4 



e 4 



-Z2 



J-i/pJ-00 

* /-'/^ /" -A(S)^2 r—-Tt yi, , , 

J-i/pJ-00 OTTp 



+ / / , 

^0 J-XjpJO 



'y{s,pz^)—vt 



-\{s)Z2 _i£p /:; 77 TTT , , , 

— e 4 -^1 + 7'(s, 02:2 a2;i as 



In the sequel, we shall need the asymptotic properties of the arc-length 
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term + 7'(s, pziY- Taylor's expansion at s = 0, z = gives 

Vi + y (0,0)2 ^i + y (0,0)2 

where ^1,^2,^3 are functions that are assumed to be smooth and bounded 
on the segment (r, ^) connecting the points (0,0) and {s,z). Hence we can 
write 



'l + 7'(s, zi^/t^sf = 1 + 0(s2 + (t - s)zl). 
In the same way, 

7(s, z) = 7(0, 0) + 7t(0, 0)s + 7'(0, 0)z + ©(s' + z'') ^vs + 0(s'' + z^), 
i.e., 

7(s, ^iVt - s) = vs + 0{s'^ + {t- s)zl). 
Let us denote for simphcity 

p Vt — s 

Again, the second integral I12 is small. Since exp(— 2;|) < 1, we have 

I/12I < C f ['^' r ^^k}L\z2\e-i{l + 0(s^ + {t-s)zl))dz2dzids 
Jo J-i/p Jo yt — s 

< C / 1] e 4, {I -\- s + [t — s)Zi) dzids 

Jq J-i/p yt - s 

■e~~ (1 + + - s)zl) dzi ds 

< C 1^ ^^{{t- s) + ^){l + s^ + {t- s))ds 
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It remains to compute In. Using the Taylor expansion we obtain 



= t ^e-i{l + 0{s^+{t-s)zl))dzids 
Jo J-i/p 47rp 

= / )== I e ^ dzi ds - ^== / e dzi ds 

Jo 47rV* - s J-oo Jo 47rV^ - s J_oo 

/■* A(s) r -d, , 

— / ^= / e 4 02:1 as 

Jo 47rV^ - s Ji/p 

+ r C'' 0{s^ + {t-s)zl)e-i dzids 

Jo 47rV* - s 

= /ill + -^112 + -^113 + hlA 

It turns out that the first integral /m gives the desired value, while the 
remaining integrals are small. To see this, we write 

lo y/iTiy/t^^^ lo V47r\/t - s ^X{0) + 0{t ), 

f'mi r"'|,.|e4<i.,d.= /'m-wt.,, = 0(e-"/'). 
Jo 47r y_oo Jo 27r 



I/112I < 

The estimate for /113 is similar and 



< C f ^^kH r{s'' + {t-s)zl)e-'^dz,ds 

Jo Vt- S J-oo 



< c 1' ^^^{s^ + t-s)ds 
vt - s 

3/2 X 



= 0{t 

Gathering the results, from the relation u{t, 0, vt) = 1/2 we have arrived 
at the identity 

Therefore, 

V = -n + 2X{Q) + 0{t), t^0 + . 
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Appendix C. Derivation of conditions in the multiphase setting 



We derive the equation and free boundary conditions for the case of mul- 
tiple phases. Since the technical aspect is similar to that of the two-phase 
case, we will give only the main ideas. 



Let us denote the phase regions by Pi. 



k, and the interface 



between regions Pi and Pj by 7jj. The symbol Ui will mean the value of the 
inner product u ■ p^, where Pi is one of the reference vectors constructed in 
Section 12.2.21 The simplest situation for three phases is depicted in figure 
C.14. Here the whole sets {ui = uj}, i,j = 1,2,3 are shown, and the parts 
of these sets corresponding to the interfaces are drawn with thicker lines. 
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Figure C.14: Conditions holding on the interfaces in the case of 3 phases (left) and con- 
figuration of reference vectors (right). 



We consider an arbitrary interface 7jj. A function 
in the above figure) can be expressed as a linear combination of reference 
vectors in the following way: 

Here = {Pi—Pj)/\Pi—Pj\ is a unit vector orthogonal to Cij = span{pi, I ^ 
and /3 = Uij where Uij := u ■ p^j. Since Pij = —Pji and Uij = —Uji, we 
will use the symbols Uij and Pij only for indices i,j satisfying i < j to avoid 
confusion. 

From the condition Ui = uj on 'jij and since Ui > ui inside Pi whenever 
i ^ I, we obtain 

/3 = 0, ai<0 on 7^,-. (C.2) 
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On the other hand, using the relation • pj 
Dirichlet functional as 



J{u) 



iVul dx 



-\/{k — 1), we compute the 

(C.3) 



In view of the separateness of conditions on ai,P in flC.2p . we can proceed 
similarly as in the two-phase case to arrive at a free boundary problem cor- 



responding to the steepest descent of J under area constraints \Pi\ 
1 = 1 k: 







= Aw 


du 








= 




dn 


dn J 




= Kj 










'dui' 








dn 




= 




7ij 



in Pi, I 
on dQ, 

hi = 1, • 



/c, i < j, 



k. 



= A, 

(C.4) 
(C.5) 

(C.6) 
(C.7) 



Here [w]^.. denotes the jump of w across %j in the normal direction and Aj/s 
are suitable functions of time. 

We assert that the above derived free boundary problem arises from the 
unconstrained gradient flow of the functional 

fc-i 



iv«<i=+E^'n 



dx. 



1=1 



This assertion is natural from a formal standpoint of the theory of Lagrange 
multipliers, since only the phase areas, multiplied by Lagrange multipliers 
A;, are added to the Dirichlet integral. However, this can be also proved by 
calculations of the first and inner variations of the functional. 

We present just an outline of the calculation of (\C.6\i for a selected inter- 
face jij, i < j < k. Equations (lC.4p and (IC.Sp are immediate and the identity 
(]C.7p is derived in an analogous way as (lC.6p . The idea is to decompose u 
again as in (IC.ip and consider the following perturbation: 



iji 
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where rj^ : R*" — )■ is a function defined by 

7]^{x) = x + eCix), 

with ( : R™ — )■ R"* a smooth function supported in a neighborhood of the 
interface •jij and a positive distance away from every other set {ui = Um}- Be- 
cause of the location of the support of (, all characteristic functions Xu pi>u p^ 
in the expression for remain unaffected by the introduced perturbation 
except for the terms Xu-p^yup^ and Xup^yu-p^- Therefore, using the reformu- 
lation of the Dirichlet integral (]C.3p . we have 

f\u)-f\u') = J ||Vm,/- |VMf/ + 

l=i,j mj^l l=i,j rriy^l 

After a long technical computation, which we omit, we obtain a condition 
holding on the interface 

/ \\\Vuij\\-2{Vu,yC)VuJ ■z/+(A,-A,)(C-z/))ci/ = 0, 

where z/ = — (V'Uij)p./|(V'Ujj)p. | is the unit outer normal to phase region Pi 
on ■jij. Since ( was arbitrary, this yields the identity ( \C.6\i with Ay = Aj — \j. 



From the above results one can deduce that the nonlinear PDE corre- 
sponding to (fT2l) in the multiphase setting will be 



k-l 



ut = Au + J2 AiP.H"-^ [dP,, (C.8) 



i=l 



where the A^'s are piecewise constant on the interfaces: 

Indeed, taking the inner product of equation ( ]C.8|) with the vector we 
use the orthogonality properties and find that 



{ui^)t = Auim + (A/ - L 
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in a neighborhood of 7^^- This means that the function uim satisfies a scalar 
equation of the type ( IT^ and thus the condition (IC.6P holds for all admissible 



In the same vein, we can use the calculation from [Appendix B] to show 
that each interface moves with normal velocity equal to minus mean curvature 
plus a space-independent term, resulting in area preservation. 
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